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1 Introduction 

The field of quantum chaos is developing very fast and there has been sub- 
stantial progress in our understanding of generic properties of eigenstates in 
classically nonintegrable and chaotic bound systems. In contrast to the the- 
oretical description of the energy spectra (and of other quantal observables) 
where we have now a rather complete understanding of the spectral statisti- 
cal universality classes and also of statistics in the transition region between 
integrability and ergodicity (i.e. going from Poisson to GOE/GUE), we are 
still far from a correspondingly complete knowledge of generic and statistical 
properties of the wavef unctions. For a few recent reviews see contributions 
in Giannoni et al (1991), Gutzwiller's book (1990), the papers in Casati et 
al (1993) and also the review on statistical properties of energy spectra by 
Robnik (1994). 

In order to understand the wavefunctions especially in the semiclassical limit 
it is intuitively very appealing to use the so-called Principle of Uniform 
Semiclassical Condensation (PUSC) of the Wigner functions (of the eigen- 
states) which is implicit in (Berry 1977a): As — we assume that the 
Wigner function of a given eigenstate uniformly (ergodically) condenses on 
the classical invariant object on which the classical motion is ergodic and 
which supports the underlying quantal state. Such an object can be e.g. an 
invariant torus, a chaotic region as a proper subset of the energy surface, or 
the entire energy surface if the system has ergodic dynamics there. 

In classically integrable systems the eigenf unctions possess a lot of ordered 
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structure globally and locally: Applying PUSC the average probability den- 
sity in the configuration space is seen to be determined by the projection 
of the corresponding quantized invariant torus onto the configuration space, 
which implies the global order. Moreover, the local structure is implied by 
the fact that the wavefunction in the semiclassical limit is locally a super- 
position of a finite number of plane waves (with the same wavenumber as 
determined by the classical momentum). 

In the opposite extreme of a classically ergodic system PUSC predicts that 
the average probability density is determined by the microcanonical Wigner 
function. Its local structure is spanned by the superposition of infinitely many 
plane waves with random phases and equal wavenumber. The random phases 
might be justified by the classical ergodicity and this assumption, originally 
due to Berry (1977b), is a good starting approximation which immediately 
predicts locally the Gaussian randomness for the probability amplitude dis- 
tribution. Berry (1977b) has also calculated the autocorrelation function of 
semiclassical chaotic (ergodic) wavefunctions which we will discuss later on 
in detail. One major surprise in this research was Heller's discovery (1984) 
of scars of unstable classical periodic orbits in classically ergodic systems. 
The scar phenomenon is of course a consequence of subtle correlations in the 
quantal phases. This has been analyzed and discussed by Bogomolny (1988) 
and Berry (1989) in the context of the Gutzwiller periodic orbit theory. The 
insufficiency of the single-periodic-orbit theory of scars has been discussed 
by Prosen and Robnik (1993a) in a study of the transition region between 
integr ability and chaos. 

In the generic case of a KAM-like system with mixed classical dynamics the 
application of PUSC is again very useful and has a great predictive power. 
Here the states can be classified as either regular (they " five" on a quantized 
invariant torus) or irregular (they "live" on a chaotic invariant region), quite 
in agreement with Percival's (1973) speculative prediction, which has been 
recently carefully re-analyzed by Prosen and Robnik (1994a). In this case 
PUSC implies asymptotic {h — > 0) statistical independence of level series 
(subsequences) associated with different regular and irregular components. 
This picture has been used by Berry and Robnik (1984) to deduce the re- 
sulting energy level statistics in such generic Hamilton systems with mixed 
classical dynamics, especially the level spacing distribution. In the recent 
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work Prosen and Robnik (1994b) have numerically confirmed the applicabil- 
ity of the Berry- Robnik theory and also explained the Brody-like behaviour 
(as discovered and described in (Prosen and Robnik 1993b)) before reaching 
the far semiclassical limit. 

2 The definition of the biUiard system and 
the numerical technique 

In the present paper we study the chaotic wavcfunctions in the 2-dim billiard 
system whose domain B (in w-plane) is defined by the complex quadratic 
conformal map of the unit disk (in z- plane) , namely 

Bx ^ {w\w ^ z + Xz"^ , \z\<l}. (1) 

as introduced by Robnik (1983,1984) and further studied by Prosen and Rob- 
nik (1993b). See also (Hayh et al 1987) and (Bruus and Stone 1994, Stone 
and Bruus 1993, 1994). Following most people in the field we shall refer 
to it as Robnik billiard. As the shape parameter A changes from to 1/2 
this system goes from the integrable case of the circular billiard continu- 
ously through a KAM-like regime to an almost ergodic regime at large A. At 
A < 1/4 the boundary is convex and therefore the Lazutkin like caustics and 
invariant tori (of boundary glancing orbits) exist. At A > 1/4 the biUiard 
was speculated (based on numerical evidence in (Robnik 1983)) to become 
ergodic, which has been disproved by Hayli et al (1987): Close to A > 1/4 
there are still some stable periodic orbits surrounded by very tiny stability 
islands. On the other hand, for A = 1/2 (the cardiod billiard) the ergodicity 
and mixing have been rigorously proved by Markarian (1993). Nevertheless, 
at large values of A, say A = 0.375 (which we study exclusively in the present 
paper) the numerical evidence does not exclude the possibility of ergodicity: 
If there are some tiny regions of stability, then they must be so small that 
they cannot be detected at large scales. 

We want to calculate and analyze the high-lying states far in the semi- 
classical limit, as high as 100,000th eigenfunction (of even parity which is 
about 200,000th when counting all states) and above, in the regime where 
the classical dynamics is almost completely ergodic (within the numerical 
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resolution of the Poincare Surface Of Section). As mentioned above, the lat- 
ter condition is satisfied at A = 0.375. However in order to reach the said 
high-lying eigenstates using the available supercomputer facilities we had 
to abandon the conformal mapping diagonalization technique developed in 
(Robnik 1984) and further employed by Prosen and Robnik (1993b). Instead 
we have implemented the Heller's method of the plane wave decomposition 
of the wavefunctions (e.g. see Heller 1991). Heller's method enables one to 
go very high in the semiclassical limit (high energies) where we can then 
calculate a few consecutive levels, whereas the diagonalization method (with 
the conformal mapping technique) has the advantage of yielding many lev- 
els from the ground state upwards. So, if one is interested in significant 
statistical analysis the latter method is superior, whilst when studying the 
individual high-lying eigenstates the former method is the better one. 

Let us spend just a few words on the technical aspects of this difficult task, 
since to the best of our knowledge many crucial ingredients have not been 
discussed in the literature so far. To solve the Schrodinger equation with 
Dirichlet boundary condition, 

A^' + E^' = 0, ^ = at the boundary, (2) 

we use the superposition of plane waves with the wavevectors of the same 
magnitude k but with different directions. The wavefunction we used for the 
even parity is 

TV 

\E'(-u, f ) = ^ CLj cos{kjuU + (pj) cos{kj^v), (3) 

where kju = kcos{9j), kjy = ksm{9j), k"^ = E the eigenenergy, N the 
number of plane waves and (pj are random phases drawn from the interval 
[0,27r), assuming uniform distribution, and 9j = 2j7T/N (i.e. the direction 
angles of the wavevectors are chosen equidistantly). The ansatz (|^) solves 
the Schrodinger equation (^ in the interior of the billiard region, so that we 
have only to satisfy the Dirichlet boundary condition. Taking the random 
phases, as we discovered, is equivalent to spreading the origins of plane waves 
all over the billiard region, and at the same time this results in reducing the 
CPU-time by almost a factor of ten. For a given k we put the wavefunc- 
tion equal to zero at a finite number M of boundary points (primary nodes) 
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and equal to 1 at an arbitrarily chosen interior point. Of course, M > N. 
This gives an inhomogeneous set of equations which can be solved by ma- 
trix inversion. Usually the matrix is very singular, thus the Singular Value 
Decomposition (SVD) method has been invoked (Press et al 1986). After ob- 
taining the coefficients aj we calculate the wavefunctions at other boundary 
points (secondary nodes). The sum of the squares of the wavefunction at 
all the secondary nodes (Heller called this sum "tension") would be ideally 
zero if /c^ is an eigenvalue. In practice it is a positive number. Therefore 
the eigenvalue problem now is to find the minimum of the "tension". In our 
numerical procedure we have looked for the zeros of the first derivative of 
the tension; namely the derivative is available analytically/explicitely from 
(^ once the amplitudes Qj have been found. In fact, since the SVD-method 
is based on finding the least square solution of the linear equations, we can 
choose M larger than N without running into the overdetermination prob- 
lem. This has been done indeed, with a typical choice M = 5N/3. It must 
be pointed out that the wavefunctions obtained in this way are not (yet) 
normalized, due to the arbitrary choice of the interior point where the value 
of the wavefunction has been arbitrarilly set equal to unity. We therefore 
explicitly normalize these wavefunctions before embarking to the analysis of 
their statistical properties. 

The accuracy of this method of course depends on the number of plane waves 
(A^) and on the number of the primary nodes (M), and we have a considerable 
freedom in choosing and M > N. In order to reach a sufficient accuracy 
the experience shows that we should take at least A^ = 3C / XdeBrogiie, and 
M = 5N/3, where C is the perimeter of the billiard and XdeBrogiie is the de 
Broglie wavelength = 2n/k. With this choice we reach the double precision 
accuracy (sixteen digits) for all levels of integrable systems like rectangular 
billiard (where the eigenenergies can be given trivially analytically) and the 
circular billiard, but also for the billiard Bx for A < 0.2. Also, the same 
choice enabled us to calculate the 100,000th even parity eigenstate and a few 
nearby eigenstates for our billiard at A = 0.375 within an accuracy of 1% of 
the mean level spacing (seven valid digits). These accuracy checks were based 
on very careful selfconsistent checks of the method and also on comparison of 
the eigenvalues with those obtained by using Robnik's diagnalization method. 

The advantage of this method is that, at one side, it is very fiexible to 
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calculate the eigenvalues, and on the other hand, it is self-checkable: The 
accuracy and the reliability can be checked by changing the interior point 
and by changing as well as M. The drawback of the method is that 
with unlucky choice of the interior point and unlucky energy step size some 
eigenstates may be — and typically are! — missed, so that the calculation 
must be repeated by using different interior points to finally collect all the 
levels. The Weyl formula (with perimeter and curvature corrections) can be 
used to detect the missing of levels (c.f. Bohigas 1991). A similar numerical 
experience has been reported in (Frisk 1990). 

3 The wavefunctions and the probabihty am- 
phtude distribution 

All the wavefunctions that we calculated and discuss here are the even parity 
eigenstates of the billiard at A = 0.375. In figure 1 we plot the even parity 
eigenfunction of energy E = 625084.5, which is about 100,010th eigenstate 
of even parity, as estimated by using the Weyl formula (with perimeter and 
curvature corrections). 



N (E) = '-±^E - + 2A)E(v/8A/(l + 2A)) - 1 ^ _ 1 

where E(x) is the complete elliptic integral of second kind (c.f. Prosen and 
Robnik 1993b). 

This is a good example of a chaotic quantum eigenstate, which exhibits 
the characteristic filamentary structure as noticed already by Heller et al 
(1987,1991), which is a consequence of the fact that in the ansatz (^ all plane 
waves (with random phases) have the same magnitude k of the wavevector. 
Also, as judged by the naked eye, the average probability density is constant 
only if the local averaging region is sufficiently large in units of de Broglie 
wavelength: Probably we need the typical size of at least several ten wave- 
lengths. The local and global average value < > is theoretically expected 
to be equal to where A = 7r(l + 2A^) is the area of B\. This is a 

direct consequence of the microcanonically uniform Wigner function for this 
eigenstate (Berry 1977, Voros 1979, Shnirelman 1979, see also Berry 1983), 
which in turn is a consequence of PUSC as explained and discussed in the 
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introduction. Indeed, the theoretical value of < > is 0.24844, whereas 
the numerical evaluation yields 0.24832 (after averaging over 1,145,294 grid 
points distributed uniformly inside the interior of the billiard region), which 
can be considered as an excellent agreement. In table 1 we compare the 
theoretical values and the numerical estimates for a number of eigenstates 
to show that quite generally this agreement is very good. There we give the 
numerical estimate for all the lowest four moments of \E'-distribution, namely 
the average mi =< \I' >, the variance m2 =< (\E' — mi)^ >, the skewness 
ma =< (\E' — mi)^ > Indl"^ and the kurtosis =< (\I' — mi)^ > /ml — 3. 
These experimental values are compared with the theoretical values of the 
Gaussian random model (see introduction) which predicts 

where again according to PUSC = 1/A = l/(vr(l + 2A^)), and ideally 
it should be equal to m2- In figure 2(a-b) we plot the numerical histogram 
for P(\l/) (merely for illustrative purposes), and — more importantly - the 
cumulative distribution /(^I') = J"^^ P{x)dx which is compared with the 
theoretical model (^). The agreement is seen to be excellent, even in figure 
2(b) where, as we shall see and discuss in section 5, in the probability density 
plot of the underlying wavefunction a scar-like feature is observed. 

In addition, we have estimated the significance levels of the cumulative 
distribution /(^E') according to the Kolmogorov-Smirnov test (Press et al 1986 
p472) with respect to the Gaussian distribution for all eigenstates listed in 
table 1. It is found to be exactly 1 within five digits. This shows again that 
indeed our results agree excellently with the theoretical prediction. 

Similar study of chaotic eigenfunctions has been published in Chirikov 
et al (1989), where even the differences between the numerics and Gaussian 
random model due to the finite dimensionality of the system have been seen. 

Our results are comparable to the findings of Aurich and Steiner (1993) 
who studied the chaotic wavefunctions of the quantum system whose classical 
counterpart is the geodesic motion on a compact surface of constant negative 
curvature, although with our numerical wavefunctions we are considerably 
farther in the semiclassical limit. So far we have not found any examples of 
scars in these high-lying states around 100,000th. (However, see section 5.) 
The conclusion is that scars are difficult to find since they "live" on smaller 
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and smaller support as fi ^ 0, or E ^ oo, and consequently asymptotically 
no longer influence the -P(^l/) distribution, as further explained in section 5. 

4 The autocorrelation function of the wave- 
functions 

The mean statistical properties of chaotic wavefunctions have been discussed, 
analyzed and described in the previous section. However the question of 
space-correlations of a wavefunction is far from trivial. The correlations ex- 
ist on different scales and their strength can vary substantially. For example, 
in figure 1 we clearly see that there is some kind of clustering on the scale of 
a few ten de Broglie wavelengths: there are regions of this size with enhanced 
probability density, and there are also regions of this size with notably de- 
pleted probability density (holes). Not every state is like that and in figure 3 
we show another chaotic even parity eigenstate with energy E = 625, 118.4, 
approximate number Ng^ = 100, 015, again for A = 0.375. Here the above 
mentioned clustering is much less pronounced and this property is well cap- 
tured in the autocorrelation function of the eigenstates as we shall see in a 
moment. 

The definition of the autocorrelation function of the probability amplitude 
of a given eigenstate is (Berry 1977b, 1983) 

C(X; q) =< vl>(q + X/2)^*{q - X/2) > / < > (6) 

where the local average denoted by < ■ ■ ■ > is taken over sufficiently large 
region around q whose size is typically many de Broglie wavelengths but still 
small compared with the geometrical size. In our case the wavefunction is 
of course real, i.e. \E'* = \E'. It should be noted that the nominator in (|^) is 
actually the Fourier transform of the Wigner function iy(q, p), 

W(q,p) = ^ /ci'Xexp(-zp-X)vI/(q-X/2)vI/(q + X/2) (7) 

where we have specialized to our real \Ef case, and also two degrees of freedom 
and h = 1. Now using the PUSC for a chaotic state following Berry (1977b) 
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we assume that the Wigner function of such a classically ergodic state is 
microcanonical, i.e. 



Wiq,p) 



6{E-H{ci,p)) 



(8) 



Jd'^qd^p6{E - H{q,p)) 



Thus substituting (H) into the inversion of (|^) and then into @ we immedi- 
ately obtain the special case of Berry's (1977b) result, namely 



where Jo is the Bessel function of zero order, k"^ is the eigenenergy and s 
is the length of X. So the autocorrelation function is isotropic, and we are 
going to check numerically the validity of this theoretical prediction. 

First we would like to check the isotropy of the autocorrelation function. 
To this end we have evaluated (^) by taking the local average on a small 
strip of 20 X 100 wavelengths situated at the center of the billiard as far as 
possible from the boundaries. The results for the wavefunctions of figure 1 
and figure 3 are shown in figures 4 and 5 correspondingly. Because of the in- 
version symmetry of the autocorrelation function w.r.t. X and the reflection 
symmetry of the wavefunctions \1/ w.r.t. v we can restrict ourselves to the 
angles within the interval [0, 7r/2], and we have chosen the values between 
and n/2 in equal steps of 7r/12, as indicated in the upper right corner of the 
figures. The autocorrelation function is obviously strongly direction depen- 
dent (please notice that the statistical noise is practically zero) and in the 
case of more uniformly chaotic wavefunction of figure 3 agrees better with 
the theoretical prediction (P) than for the less chaotic eigenstate of figure 1. 
We believe that the semiclassical periodic orbit theory (see e.g. Casati et al 
1993, Tel and Ott 1993 and references therein) could explain the deviations 
from the isotropy. Our results agree qualitatively with (Aurich and Steiner 
1993) although we are considerably higher in the semiclassical limit (by a 
factor of 10 or so), and also with somewhat old results in (McDonald and 
Kaufman 1988, Shapiro and Goelman 1984, Shapiro et al 1988). 

It is interesting that after averaging over many directions we get a consider- 
able agreement with (|^). This is shown in figure 6(a-d) where we vary the 
size of the averaging disk and also the number of the directions over which 



C(X;q) = Jo(M 



(9) 
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the average is taken. These plots are for the eigenfunction shown in figure 
1. Both effects are clearly visible, namely the increasingly better agreement 
with as we increase the radius of the averaging disk and/or as we increase 
the number of directions. The same aspects are shown in figure 7(a-d) for 
the more uniformly chaotic state of figure 3. By comparing the figures 6(a) 
and 7(a) we see that in the latter plot the agreement with theory is better. 

5 Scar-like features in wavefunctions 

As we know since Heller's (1984) discovery of scars (of unstable classical peri- 
odic orbits) in chaotic quantum eigenf unctions of classically ergodic systems 
we do expect such scars to exist in all chaotic systems, but according to 
the single-periodic-orbit theory (Bogomolny 1988, Berry 1989) the scar sup- 
porting region should shrink as \fTi as ^ — >• 0, whilst the probability density 
contrast remains fixed since it is predicted to be ^-independent (Heller 1984). 
The many-orbits theory (Robnik 1989) would speculatively predict the linear 
scaling of the with ^ as a consequence of the interference effects. 

Some phenomenological material on this topic has been recently published 
in (Prosen and Robnik 1993a). 

In this short section we would like to draw attention to an interesting scar-like 
feature seen in figure 3: Near the boundary, about ten de Broglie wavelengths 
away, there is a thin scar-like feature which has no simple explanation, be- 
cause due to the nonconvexity of the billiard boundary there are no Lazutkin- 
like caustics and invariant tori and also no such glancing periodic orbits. The 
only classical object that might be relevant for this feature is possibly the 
glancing orbit which survives many bounces while going round the boundary 
until reaching the non-convexity region and fiying away, becoming completely 
chaotic afterwards. We have observed a few similar features in quite a few 
other eigenfunctions, but we cannot offer any definite theoretical explanation 
so far. However, the formalism offfered and discussed in (Miiller et al 1993) 
might just be right to quantitatively describe the role of such orbits which 
are recurrent in configuration space but not periodic. 
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6 Discussion and conclusions 



In this paper we have numerically calculated the high-lying chaotic states 
in the Robnik billiard as high as 200,000th eigenstate and investigated their 
semiclassical morphology and their statistical properties. To achieve this we 
have implemented and adapted Heller's method of plane wave decomposition 
which has been further developed and its accuracy carefully checked. Similar 
to other workers (e.g. Aurich and Steiner 1993) we reach the following con- 
clusions. In such high-lying eigenstates the scars are hardly detectable since 
so far we have not found any of them. The average probability density is glob- 
ally in excellent agreement with the theoretical semiclassical (and classical!) 
prediction. The Gaussian random model for the local statistical properties 
of the wavef unctions is generally excellent, in spite of the characteristic fila- 
mentary structure and the relevant clustering of probability density on the 
scales of a few ten de Broglie wavelengths. This has been found by comparing 
the theoretical and the numerical distributions and also by the comparison of 
the lowest four moments and the evaluation of the Kolmogorov-Smirnov test. 
The autocorrelation function nicely captures the clustering property, is found 
to be strongly direction dependent in contradistinction with Berry's (1977b) 
isotropic prediction, but after averaging over many directions the agreement 
with Berry's theory is recovered. Finally we should mention that in some 
of the eigenstates we discovered scar-like features resembling the whispering 
gallery modes for which we do not have a proper theoretical explanation. 

Our current and future work deals with the systematic search for the scars 
and the analysis of their geometry and scaling properties with h. On theo- 
retical side the present paper stimulates further work on the scar theory for 
which we expect improvement when many-orbits theory will be set up follow- 
ing the suggestions in (Robnik 1989, Prosen and Robnik 1993a). Moreover, 
we believe that the application of the Gutzwiller's (one-) periodic orbit the- 
ory could explain in detail the anisotropics of the autocorrelation function. 
Our work also shows that there is still much interesting structure in the range 
of a few ten de Broglie wavelengths in chaotic wavefunctions which calls for 
a more refined statistical description. 
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Table 



Table 1: The average, variance, skewness and kurtosis of a few eigenstates 
nearby 100,000th eigenstate of even parity in comparison with theoretical 
values. Nev{E) is given by the Weyl formula (^. The significance levels 
of the Kolmogorov-Smirnov test for all eigenstates listed in this table are 
exactly one within 5 digits. 



E 


Ne.{E) 


Average 


Variance 


Skewness 


Kurtosis 


625040.6 


100003 


0.00002 


0.24^ 


^28 


-0.00211 


0.11607 


625058.4 


100006 


0.00001 


0.24^ 


^36 


0.00441 


0.06576 


625084.5 


100010 


0.00001 


0.24^ 


^32 


-0.00160 


0.06809 


625099.5 


100012 


-0.00001 


0.24^ 


^34 


-0.00179 


0.03638 


625118.4 


100015 


0.00001 


0.24^ 


^38 


0.00460 


-0.01415 


625161.9 


100022 


0.00003 


0.24^ 


^37 


0.00347 


0.02424 


625172.8 


100024 


0.00007 


0.24^ 


^39 


-0.00248 


0.04836 


625182.1 


100025 


0.00000 


0.24^ 


^54 


-0.00135 


0.03624 


Gaussian 




0.0 


0.24^ 


M4 


0.0 


0.0 
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Figure captions 



Figure 1: The probability density plot for the even parity eigenstate with 
E = 625, 084.5, for A = 0.375, and with the estimated sequential number 
using the Weyl formula equal to 100,010. The contours are plotted at ten 
equally spaced steps between zero and the maximum value. In this geometry 
the unit length is 126 de Broglie wavelengths. 

Figure 2: The probability distribution function and cumulative distribution 
function of eigenstates E = 625, 084.5 (approximately 100,010th even parity 
state) (a), and E = 625, 118.4 (approximately 100,015th even parity state) 
(b), in comparison with the Gaussian random model (^. In the top dia- 
grams we show the histograms compared with theoretical curve (^, and in 
the lower diagrams we show the cumulative amplitude distribution function 
/(\E'). Three small boxed regions are displayed in the corresponding magni- 
fied windows. Here the difference between the theoretical and the numerical 
curves is hardly visible since the agreement is so good. 

Figure 3: The probability density plot for the even parity eigenstate with 
E = 625, 118.4, for A = 0.375, and the estimated sequential number using the 
Weyl formula (^ equal to 100,015. The contours are plotted at ten equally 
spaced steps between zero and the maximum value. In this geometry the 
unit length is 126 de Broglie wavelengths. 

Figure 4: The autocorrelation function C(X, q) of the eigenstate in figure 
1. C(X, q) is plotted against ks for seven different angles of X w.r.t the 
abscissa; here s = |X|. The angle and the averaging strip are indicated in 
the upper right corner of the figure. From (a) to (g) the angle goes from 
to 7r/2 with the increment of 7r/12. The dashed curve is the theoretical 
prediction (^, namely Jo{ks) , whilst the full curve denotes the numerical 
result. The local average has been taken on a strip of 20 x 100 de Broglie 
wavelengths. For a fixed ks about 150,000 grid points inside the strip have 
been used to calculate C. The reference point q is fixed at (0.4, 0.0) which 
is probably sufficiently far away from the billiard boundary. 
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Figure 5: The same as figure 4 but for tlie eigenstate of figure 3. 

Figure 6: Tlie autocorrelation function after averaging over many directions 
X of tlie eigenstate in figure 1. C(X, q) is plotted against ks for three 
different averaging disks and different number of directions, where s — |X|. 
The averaging disk in (a), (b), (c) and (d) has the diameter of 200, 100, 50 
and 100 de Broglic wavelengths, respectively. In (a), (b) and (c) the number 
of directions is 200, whilst in (d) it is 400. 

Figure 7: The same as in figure 6 but for the eigenstate in figure 3. 
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